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Abstract 

It has been recently pointed out that dynamical systems depending on future values of the unknowns may 
be useful in different areas of knowledge. We explore in this context the extension of the concept of order 
reduction that has been useful with singular and delay differential equations in electrodynamics and general 
relativity. We discuss some general properties of order reductions in this new context and explore a method 
of successive approximations, which among other things is used to check and improve the "extrapolate 
prediction" and "fixed rate prediction" methods. 
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I. INTRODUCTION 



T. Ohira has recently proposed and analyzed a formalism and concrete examples of dynamical 
systems governed by predictions of future states, which he calls "predictive dynamical systems" 
111]. Since initial conditions are not sufficient to solve this kind of dynamical system (and make 
sure the solution is unique), Ohira proposes two methods of predicting the future values of the 
unknowns necessary to find (numerically) the solution to the system: "fixed rate prediction" and 
"extrapolate prediction." Since both methods can at most provide some approximate solution of 
the dynamical system, it may be interesting to explore other methods that could eventually improve 
the quality of the approximation. 

In classical electrodynamics and general relativity one finds singular differential equations and 
delay differential equations for which the usual physical initial conditions are not enough to com- 
pute the solution. In this context the idea of order reduction has been useful Isl IsiH. Q. S S 

We extend the concept of order reduction to predictive dynamical systems in section lUl as well 
as a method of successive approximations to compute the solution. Since we lack a rigorous theory 
of this method, a simple but illustrative example is analyzed in section Uni Numerical results are 
discussed in section |lVl for the same dynamical systems discussed in jl|. In section |V| we briefly 
explore higher order reductions and compare again our numerical results with those of Ohira's. 



II. ORDER REDUCTIONS OF PREDICTIVE DYNAMICAL SYSTEMS 

Although our definition and results can be readily extended to continuous dynamical systems, 
for simplicity we will consider only discrete dynamical systems in which a physical quantity x is 
defined only for integer values ?7, = 0, 1, 2,. ..of "time" according to a law in the form 

Xn+l = M Xn+p) , (1) 

with some advance p = 1, 2, ... It is obvious that an initial condition xq for, say, tt, = is not 
enough to predict the future. Even if one could solve ([T]) for x„+p, the resulting dynamical system 

^n+p = N i^m ^n+l) (2) 

would require specifying p initial conditions: xq, xi, . . . , Xp_i. In consequence more assumptions 
are necessary to solve (HJi. In Ohira's "fixed rate prediction" method 111, one replaces Xn+p on the 
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right hand side of (HJ) by + p {xn — Xn^i), which still needs some additional initial condition 
for n = —1 or another suitable assumption. In the "extrapolate prediction" method one would 
substitute for Xn+p the value obtained by applying p times to Xn the map ([T]) with p = 0. It is clear 
that, in general, both methods will provide at most approximations to a solution of the original 
dynamical system. 

The idea behind order reductions is that dU) is not the true evolution equation but only a neces- 
sary condition every solution to the actual (unknown) dynamical system must satisfy. If the true 
dynamical system is a deterministic one in the form 

Xn+l = F {Xn) , (3) 

knowledge of the future is only necessary because our incomplete theory did not led us to © 
but only to a less restrictive condition in the form O. Since the latter must be satisfied by every 
solution to (jSj), we have the following condition for the unknown F: 

p times 

F{x) = M {x, FP{x)) , FP =FoFo---oF . (4) 

Of course, in general, one cannot solve © for F (this is the reason it is unknown) but one can 
try finding good approximations by different methods. In some cases there is a small parameter 
in the problem, so that the natural way would be to try Taylor expansions with respect to that 
parameter. But we are here going to explore a general method of successive approximations which 
have proved useful with singular al delay differential equations 

We will construct a succession of approximations Fq, Fi, . . . defined by 

F^+i(x) = M(x,F^(x)) (5) 

along with some suitable initial Fo(x). It is clear that if the succession is convergent, its limit 
F(x) = limm-^oo Fm{x) is a solution of (|4|). One obvious choice for the initial condition is 

Fo{x) = M {x,x) , (6) 

in which case Fi is the value Xn+i obtained by means of Ohira's "extrapolate prediction;" but, 
although the limit F (x„) will be unattainable in practice, the approximation can be improved by 
computing successive Fm (xn) until \Fm+i (xn) — F^ {xn) \ is below some tolerance value. How- 
ever, we will see later that different initial conditions may change dramatically the convergence 
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rate (moreover, the may lead to different order reductions), so that in practice some additional 
criterion must be used (for instance, in electrodynamics one can use the limit in which the charge 
vanishes to select the right order reduction). 

The problem of the existence of the limit F is here posed in too general grounds to have an 
answer. Instead of that we will consider an artificial but illustrative problem. 



III. A LINEAR EXAMPLE 



Let us first consider the discrete dynamical system 

+ bxn+i, (07^0,67^0,1). (7) 

Of course, this can be written as 

a 

Xn+l = OLXn, « = " (8) 

1 — 

but let us pretend we do not know that and want to solve © by the method of successive ap- 
proximations. It is easy to see that for Fq{x) = (a + b)x, which corresponds to Q, or for any 
Fq{x) = OqX with constant ao, we have 

Fmix) = amX (9) 

with 

Om+i = a + bam- (10) 

Since (a — ftm+i) = b{a — am), whatever ao is, the recurrence (fTOl) will converge (to a) if and 
only if |6| < 1. In consequence, in this example the method of successive approximations will 
converge (to the right dynamical system) when |6| < 1 and diverge for |6| > 1. One cannot 
expect the method to be convergent always, but the example suggests that (as is often the case in 
electrodynamics it may work if some parameter in the theory is small enough. 

The following example is 

Xn+i = aXn + bxn+2, {a ^ , b ^ , -1 - a) . (11) 

Also in this case we can solve for x„+2 to obtain a two-point recurrence which needs two initial 
conditions (say x_i and xq). Instead, we seek an order reduction (jSj) which only requires one 
initial condition and must satisfy 

F{x) = ax + bF{F{x)). (12) 



For Aab < 1 this functional equation has, at least, the following two linear solutions: 

r( \ _ 1±V1-Aab 

r [X) = ax, a = a± = ; . (13) 

26 

Starting from any Fo{x) = aox with constant we get again Q with 

am+i = a + ba'^. (14) 

But this quadratic map is just the logistic map whose properties have been explored in depth in 
chaos theory [12]. For this reason it is easy to prove that a-m will converge to a_ for any parameter 
values such that —3 < 4a6 < 1 provided the initial condition is choosed so that 



l«o| < . (15) 

This is the case for ao = a + b — which correspond to Q — for \a + b\ small enough. For other 
initial conditions or parameters may go to infinity, approach a cycle of any period or change 
chaotically. Again we see that the method could work for small parameter values, but also that 
it could never converge to the right solution iF{x) = a^x, for instance), in which case other 
methods should be tried (maybe an appropriate series expansion, or a numerical method to solve 
O for Xn+p). 



IV. NUMERICAL RESULTS 



Successive approximations to the order reduction can be numerically computed in any program- 
ming language. For instance, the Mathematica code 1 13] in Table |I| will compute and display Xn 
(for n = 0, 1, . . . , 10) by using the second approximation F2, in the case of the "sigmoid function" 
discussed in reference yj]: 

M(x, y) = {l- + — ^ - 1. (16) 

We have use that code for Figure [T} where the values xq, xi, . . . , xio obtained with Fi, F2, F4 
and F5 are displayed for fi = 0.5, (3 = 0.8, p = 5 and initial guess ®. The dots in the upper 
polygonal have been computed with Fi and, thus, are the same obtained by Ohira's "extrapolate 
prediction." We can see there is room for improvement, for the values with F2 are rather smaller, 
while those obtained with F4 and F5 are indistinguishable in the figure, proving they are very near 
those one would obtain with the limit F. We can see in Figure |2l the importance of a good guess 
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for Fq: selecting Fo{x) = M (x, xq) leads to a much slower convergence and even Fg is not a good 
approximation. 

In Figure|3]one can see that convergence is faster for the "Mackey-Glass function" of reference 

M{x,y) = il-fi)x + -^, (17) 

with /i = 0.5, P = 0.8, s = 10, p = 5 and initial guess Q: solutions with F2 and F3 are already 
very close. 

One can also have the program compute at each step Xn successive approximations F^ 
until the difference between two consecutive approximations is below some maximum relative 
error, which is called tol in the code in Tablellllfor the "Mackey-Glass function" of reference iQ]. 
For more complex calculations this code can (must) be improved in many ways, including a better 
storage management (here every computed value is stored) and using a compiled programming 
language. 

V. HIGHER ORDER REDUCTIONS 

To keep things simple we have reduced dl]) to the first-order dynamical system (jSj), which only 
needs xq to identify each solution. In some cases we might have theoretical reasons to think that 
the true dynamical system is of second order, 

Xn+l = G(Xn,Xn-l) , (18) 

with 

Gix,y) = M{x,G^P\x,y)), (19) 

and 

G'-'\x,y)^x, G^'\x,y)^G{x,y), G^^+'\x , y) ^ G {G^^'' {x , y) , G^^-'\x , y)) , (20) 

so that x_i and xq must be specified. Notice that in the corresponding scheme of successive 
approximations, 

Gm+i{x,y) = M(x,G^^\x,y)), (21) 
one could use Ohira's "fixed rate prediction" [1] to provide the following starting guess: 

Go{x,y) = M{x,x + p{x-y)). (22) 



We have used the code in Table HiH to compute the results in Figure |4l where the results for G2, 
Gi, Ge and Gg are displayed for fi = 0.5, P = 0.8, p = 5, xq = x_i = 0.5 and initial guess 
(E^ . Wee see that the successive approximations converge slowly to the same solution displayed 
in Figure[T] in particular this means that in this example © is also an order reduction of (fTHl . 

VI. FINAL COMMENTS 

We have extended the concept of order reduction to predictive dynamical systems and discussed 
some examples in which it can be used to construct good approximations to exact solutions of those 
systems. In particular, we have shown that a method of successive approximations may be used to 
check and improve the accuracy of Ohira's extrapolate prediction fv\. We are not claiming that the 
method will work always, but that, as happens with singular and delay differential equations, there 
may be interesting cases in which it can be used to construct solutions to predictive dynamical 
systems. lin other cases one must have to resort to other approximation scheme, such as series 
expansions, backward integration, shooting methods, (or a root finding routine to solve for Xn+p 
at each step), etc. 

To keep things simple we have only considered discrete dynamical systems; but the con- 
cepts explored here can be extended in an obvious way both to reductions of higher order and 
to differential-difference equations of advanced type (with the meaning defined in reference ll^). 
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Clear [F] 

F[m_, x_] := F[m, x] = 

M[x,Nest[F[m-l,#]&,x,p] ] 

F[0,x_] :=F[0,x] =M[x,x] 

M[x_,y_] := (l-0.5)x+2/ (l+Exp[-0.8 y] )-l 

p = 5; 

xO = 0.5; 

m = 2; 

ListPlot [NestList [F [m, #] &, xO, 10] ] ; 



(* Forget previous calculation *) 

(* Recurrence *) 

(* Initial guess *) 

(* Map *) 

(* Advance *) 

(* Initial condition *) 

(* Approximation *) 

(* Plot F_m(x_n) *) 



TABLE I: Mathematica program to compute and display (a;„). 
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Clear [F] (* Forget previous calculation *) 

F :: "iterations" = "Too many iterations."; 
F[x_] := F[x] = 

Module [ {m} , (* Successive approximations *) 

For [m = 1 , 

Abs [ (F [m, x] -F [m-1, x] ) / (F [m, x] +0 . 001) ] > tol, 
m++. 

If [m > mmax. Message [F :: "iterations" ] ; Break [] ] 



F[m,x] 



F[m_, x_] := F[m, x] = 

M [x. Nest [F [m-1, #]&,x,p] ] 

F [0,x_] := F [0,x] = M[x,x] 
M[x_,y_] := (l-0.5)x+0.8y/ (1+y'lO) 
p = 8; 
xO = 0.5; 
tol = 10^-5; 
mmax = 10; 

ListPlot [NestList [F, xO, 10] ] ; 



(* Recurrence *) 

(* Initial guess *) 

(* Map *) 

( * Advance * ) 

(* Initial condition *) 

(* Maximum relative error *) 

(* Maximum value of m *) 

(* Plot F{x_n) *) 



TABLE n: Mathematica program to compute and display an approximation to F 
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Clear [G,x] (* Forget previous calculation *) 

G [ 0, m_, x_, y_] := x (* Recurrence *) 

G[p_,m_, x_, Y_] := G[p,m,x,y] = G [ 1 , m, G [p-1 , m, x, y] , G [p-2 , m, x, y ] ] 
G [l,m_, x_, y_] := G[l,m,x,y] = M [ x, G [p, m-l , x, y ] ] 



(* Initial guess *) 
(* Map *) 

( * Advance * ) 

(* Initial conditions *) 



G[l,0,x_,y_] := G[l,0,x,y] = M [x, x+p (x-y ) ] 
M[x_,y_] := (l-0.5)x+2/ (l+Exp[-0.8 y] )-l 
p = 5; 

x[m_, -1] := 0.5; 
x[m_, 0] = 0.5; 

X [m_, n_] := x[m, n] = G [ 1 , m, x [m, n-1 ] , x [m, n-2 ] ] (* Order reduction * 
m = 1 ; 

ListPlot [Table [x [m, n] , {n, 0, 50 } ] ] ; 



(* Approximation *) 
(* Plot solution *) 



TABLE III: Mathematica program to compute and display Gm (x^, Xn-i)- 
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